# 9 January 2021
# Baris Ari
# Replication of Figure 6

cat("\014")  
remove(list = ls())
library(stats)
library(dplyr)
library(ROCR)
library(survival)
library(gdata)
library(texreg)
library(MASS)


# Preliminaries ####

##Color scheme####
percent <- 60
fark <- 60
col.neg <- rgb(red = col2rgb("gray40")[1], 
               green = col2rgb("gray40")[1],
               blue = col2rgb("gray40")[1], max = 255,
               alpha = (100-percent)*255/100)
col.both <- rgb(red = col2rgb("gray60")[1] - fark*2, 
                green = col2rgb("gray60")[1] - fark,
                blue = col2rgb("gray60")[1], max = 255,
                alpha = (100-percent)*255/100)

col.ref <- rgb(red = col2rgb("gray80")[1] - fark, 
               green = col2rgb("gray80")[1] - fark/2,
               blue = col2rgb("gray80")[1], max = 255,
               alpha = (100-percent)*255/100)

col.white <- rgb(red = col2rgb("white")[1], 
                 green = col2rgb("white")[1],
                 blue = col2rgb("white")[1], max = 255,
                 alpha = (100-0)*255/100)
col2rgb("indianred")

ct <- function(x){x*255/100}


mycol <- "lightsalmon1"
col2rgb(mycol)
mycolline <- "lightcoral"

# Preliminaries  continued
mygrey  <- rgb(t(col2rgb("grey80")), max = 255, alpha = 125)
mygrey1 <- rgb(t(col2rgb("grey60")), max = 255, alpha = 125)
mygrey2 <- rgb(t(col2rgb("grey70")), max = 255, alpha = 125)
mygrey3 <- rgb(t(col2rgb("grey80")), max = 255, alpha = 125)
mygrey4 <- rgb(t(col2rgb("grey90")), max = 255, alpha = 125)

# mode finder function
mode.proper <- function(x){
  x <- na.omit(x)
  md<- as.numeric(names(table(x))[table(x)==max(table(x))])
  return(md)}


# Data ####
d <- read.csv("data.csv")

d <- d %>% dplyr::select(SideA, SideB, DyadId, ConflictId, GWNoA, Year, 
                         Negotiation, demref, 
                         polyarchy, duration_dy_ln, imp.gdppcln, imp.popln,
                         Intensity, 
                         Region, RebStrBin.imp, Ethnic.incomp, Incompatibility, 
                         medprev2y,
                         ExChange, oda.pc.ln, UNPKO.total.ln, BdBest_ln, tsln, tsln2, tsln3,
                         demref_poly_xp2_ip2_lag,
                         demref_lexical_autlib_dec5_lag, demref_lexical_demtrans_dec5_lag)

factors <- c("RebStrBin.imp", "Region", "Ethnic.incomp", "ExChange", "Intensity", "medprev2y")
table(factors %in% names(d))

for (f in factors){
  d[,f] <- as.factor(d[,f])
}


#variable renamed
d$exit <- as.numeric(as.character(d$ExChange))

# dyad duration variable renamed
d$dydurln <- d$duration_dy_ln

## Model formula:
d3 <- d %>% dplyr::select(DyadId, Year, GWNoA, Negotiation, polyarchy, demref, dydurln , 
                          imp.gdppcln, imp.popln,
                          Intensity,
                          tsln, tsln2, tsln3, RebStrBin.imp, Ethnic.incomp, exit, Region,
                          medprev2y)
d3 <- na.omit(d3)

m3    <- as.formula("Negotiation ~ polyarchy + demref + dydurln + imp.gdppcln + imp.popln +
                            Intensity + RebStrBin.imp  + Ethnic.incomp + exit + medprev2y + 
                            Region + tsln + tsln2 + tsln3
                            ")

## without Indonesia estimate
indCode  <- unique(d$GWNoA[d$SideA == "Indonesia"])
d3.WOind <- d3[d3$GWNoA != indCode, ]
m3.est <-  glm(formula = m3, data = d3.WOind, family = binomial(link = "logit"))


d3.GAM <- d[d$SideB == "GAM", ]
d3.GAM <- d3.GAM %>% dplyr::select(Year, Negotiation, polyarchy , demref , dydurln , 
                                   imp.gdppcln , imp.popln ,
                   Intensity , RebStrBin.imp , Ethnic.incomp , exit , medprev2y , 
                   Region , tsln , tsln2 , tsln3)

GAM.alternative <- d3.GAM
GAM.alternative$polyarchy <- 0.20

X <- cbind(Intercept = 1,
            polyarchy = d3.GAM$polyarchy,
            demref = d3.GAM$demref,
            dydurln = d3.GAM$dydurln,
            imp.gdppcln = d3.GAM$imp.gdppcln,
            imp.popln = d3.GAM$imp.popln,
            Intensity1 = ifelse(d3.GAM$Intensity == 1,1,0),
            Intensity2 = ifelse(d3.GAM$Intensity == 2,1,0),
            RebStrBin.imp1 = d3.GAM$RebStrBin,
            Ethnic.incomp1 = d3.GAM$Ethnic.incomp,
            exit1 = d3.GAM$exit,
            medprev2y1 = d3.GAM$medprev2y,
            Region2 = 0,
            Region3 = 1,
            Region4 = 0,
            Region5 = 0,
            tsln = d3.GAM$tsln,
            tsln2 = d3.GAM$tsln2,
            tsln3 = d3.GAM$tsln3
)

Xa <- cbind(Intercept = 1,
            polyarchy = 0.20,
            demref = 0,
            dydurln = d3.GAM$dydurln,
            imp.gdppcln = d3.GAM$imp.gdppcln,
            imp.popln = d3.GAM$imp.popln,
            Intensity1 = ifelse(d3.GAM$Intensity == 1,1,0),
            Intensity2 = ifelse(d3.GAM$Intensity == 2,1,0),
            RebStrBin.imp1 = d3.GAM$RebStrBin,
            Ethnic.incomp1 = d3.GAM$Ethnic.incomp,
            exit1 = d3.GAM$exit,
            medprev2y1 = d3.GAM$medprev2y,
            Region2 = 0,
            Region3 = 1,
            Region4 = 0,
            Region5 = 0,
            tsln = d3.GAM$tsln,
            tsln2 = d3.GAM$tsln2,
            tsln3 = d3.GAM$tsln3
)

length(colnames(model.matrix(m3.est)))
length(m3.est$coefficients)

ystar1 <- X %*% m3.est$coefficients
ystar2 <- Xa %*% m3.est$coefficients

yhat1 <- exp(ystar1) / (exp(ystar1) + 1)
yhat2 <- exp(ystar2) / (exp(ystar2) + 1)

p <- data.frame(Year = d3.GAM$Year, real = yhat1, alt = yhat2)
plot(p$Year, p$real, xlim = range(p$Year), ylim = c(0,1), t = "l")
lines(p$Year, p$alt, col = "red")

# Simulate  ####
coef.draws <- mvrnorm(n = 10000, mu = m3.est$coefficients, Sigma = vcov(m3.est))

ystar1 <- t(coef.draws %*% t(as.matrix(X)))
yhat1  <- exp(ystar1) / (exp(ystar1) + 1)

ystar2 <- t(coef.draws %*% t(as.matrix(Xa)))
yhat2 <- exp(ystar2) / (exp(ystar2) + 1)

real.pe <- apply(yhat1, 1, quantile, probs =  0.50)
real.lb <- apply(yhat1, 1, quantile, probs =  0.025)
real.ub <- apply(yhat1, 1, quantile, probs =  0.975)

alt.pe <- apply(yhat2, 1, quantile, probs =  0.50)
alt.lb <- apply(yhat2, 1, quantile, probs =  0.025)
alt.ub <- apply(yhat2, 1, quantile, probs =  0.975)


plot(p$Year, p$real, xlim = range(p$Year), ylim = c(0,1), t = "l")
lines(p$Year, real.pe, lwd = 2)
lines(p$Year, real.lb, lwd = 2, lty =2)
lines(p$Year, real.ub, lwd = 2, lty =2)


lines(p$Year, alt.pe, lwd = 2, col = "red")
lines(p$Year, alt.lb, lwd = 2, lty =2, col = "red")
lines(p$Year, alt.ub, lwd = 2, lty =2, col = "red")

d.GAM<- d3.GAM

d.GAM$yhat.sim1.pe.GAM <- real.pe

# background colors
d.GAM$colors <- "white"
d.GAM$colors[d.GAM$Negotiation == 1 & d.GAM$demref == 1] <- col.both
d.GAM$colors[d.GAM$Negotiation == 1 & d.GAM$demref == 0] <- col.neg
d.GAM$colors[d.GAM$Negotiation == 0 & d.GAM$demref == 1] <- col.ref


m3.red  <- as.formula("Negotiation ~ polyarchy + dydurln + imp.gdppcln + imp.popln +
                            Intensity + RebStrBin.imp  + Ethnic.incomp + exit + medprev2y +
                      Region + tsln + tsln2 + tsln3
                            ")
m3est.r <-  glm(formula = m3.red, data = d3.WOind, family = binomial(link = "logit"))

Xr <- cbind(Intercept = 1,
           polyarchy = d3.GAM$polyarchy,
           #demref = d3.GAM$demref,
           dydurln = d3.GAM$dydurln,
           imp.gdppcln = d3.GAM$imp.gdppcln,
           imp.popln = d3.GAM$imp.popln,
           Intensity1 = ifelse(d3.GAM$Intensity == 1,1,0),
           Intensity2 = ifelse(d3.GAM$Intensity == 2,1,0),
           RebStrBin.imp1 = d3.GAM$RebStrBin,
           Ethnic.incomp1 = d3.GAM$Ethnic.incomp,
           exit1 = d3.GAM$exit,
           medprev2y1 = d3.GAM$medprev2y,
           Region2 = 0,
           Region3 = 1,
           Region4 = 0,
           Region5 = 0,
           tsln = d3.GAM$tsln,
           tsln2 = d3.GAM$tsln2,
           tsln3 = d3.GAM$tsln3
)
ystar.r <- Xr %*% m3est.r$coefficients
yhat.r <- exp(ystar.r) / (exp(ystar.r) + 1)

#Reported in main text
real.pe[d.GAM$Year == 1998] 
real.pe[d.GAM$Year == 1999] 
real.pe[d.GAM$Year == 2000] 

yhat.r[d.GAM$Year == 1998] 
yhat.r[d.GAM$Year == 1999] 
alt.pe[d.GAM$Year == 2000] 

round(real.pe[d.GAM$Year == 1999],2)
round(real.pe[d.GAM$Year == 2000],2)  - round(alt.pe[d.GAM$Year == 2000],2)

real.pe[d.GAM$Year == 2005]  - alt.pe[d.GAM$Year == 2005] 


# Graph  ####
#add 2005.8 to stand for the last year #
d.GAM <- rbind(d.GAM, d.GAM[nrow(d.GAM),])
d.GAM$Year[nrow(d.GAM)] <- 2005 + 8/12
real.pe <- c(real.pe , real.pe[length(real.pe)])
alt.pe <- c(alt.pe , alt.pe[length(alt.pe)])
yhat.r <- c(yhat.r , yhat.r[length(yhat.r)])

dev.off()
png("indonesia_outofsample.png", 
    units = "cm", width = 24, height = 12, res = 300)
par(mar = c(2,4,2,1), mfrow = c(1,1), mgp = c(3,0.75,0))

plot(d.GAM$Year, 
     real.pe, t = "n", 
     ylim = c(0, 1), xlim = c(1990, 2006),
     col = "black", bty = "n",
     axes = F, 
     main = c(paste("Out-of-sample predictions for Indonesia vs GAM")),
     xlab = "Years",
     ylab = "Pr(Negotiation)")
for (t in 1990:2005){
  polygon(x = c(t, t + 1, t +1, t),
          y = c(-0.5, -0.5,1,1),
          border = NA,
          col = d.GAM$colors[d.GAM$Year == t])
}
polygon(x = c(d.GAM$Year[nrow(d.GAM)], 2006, 2006, d.GAM$Year[nrow(d.GAM)]),
        y = c(-0.5, -0.5,1,1), border = NA, col = "white" )

### lines ####
lines(d.GAM$Year, real.pe, t = "l", ylim = c(0, 0.8), col = "black", lwd = 4, lty = 1)
lines(d.GAM$Year, alt.pe, t = "l", ylim = c(0, 0.8), col = "black", lwd = 2, lty = 2)
lines(d.GAM$Year, yhat.r, t = "l", ylim = c(0, 0.8), col = mycolline, lwd = 3, lty = 3)

abline( v = 2005 + 8/12, lty  = 4)
text(2006, 0.50, "Peace Agreement", srt=90, cex = 1)

# lines(d.GAM$Year[d.GAM$Year >= 1997], alt.lb[d.GAM$Year >= 1997], 
#       t = "l", ylim = c(0, 0.8), col = "black", lwd = 1, lty = 3)
# lines(d.GAM$Year[d.GAM$Year >= 1997], alt.ub[d.GAM$Year >= 1997], 
#       t = "l", ylim = c(0, 0.8), col = "black", lwd = 1, lty = 3)

#axis 1 
axis(1, 1990:2006, tck = 0.01)
axis(1, 1990:2006, labels = NA, tck = -0.01)

axis(2, las = 2, at = seq(0,1,0.2), tck = -0.02)
axis(2, at = seq(0,1,0.1), tck = -0.01, labels = NA)

## legend
polygon(x = c(1990,1991,1991,1990), y = c(0.95, 0.95, 1, 1), border = NA, col = col.neg)
polygon(x = c(1990,1991,1991,1990), y = c(0.89, 0.89, 0.94, 0.94), border = NA, col = col.ref)

text("Negotiation", adj = c(0, 0.5), x= 1991.2, y=0.975, cex = 0.8)
text("Democratization", adj= c(0, 0.5), x= 1991.2, y=0.92, cex = 0.8)
#polygon(x = c(1989,1990,1990,1989), y = c(0.88, 0.88, 0.83, 0.83), border = NA, col = mycol)
segments(x0 = 1990, x1 = 1990.95, y0= 0.855, y1 = 0.855, col = "black", lwd = 3, lty = 1)
text("Full model (dem. reform variable included)", adj = c(0, 0.5), x= 1991.2, y=0.855, cex = 0.8)
text("Reduced model (dem. reform variable excluded)", adj = c(0, 0.5), x= 1991.2, y=0.785, cex = 0.8)
text("Full model counterfactual (no democratization)", adj = c(0, 0.5), x= 1991.2, y=0.715, cex = 0.8)
segments(x0 = 1990, x1 = 1991, y0= 0.785, y1 = 0.785, col = mycolline, lwd = 2, lty = 3)
segments(x0 = 1990, x1 = 1991, y0= 0.715, y1 = 0.715, col = "black", lwd = 2, lty = 2)
dev.off()


